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Phase change material To quantify their technical and economic feasibility for building’s applications, computational models 


Building enclosure 


i of TES that can be integrated into whole building energy simulations are highly demanded. This paper 
Numerical models 


reviews the different modeling methods generally used for PCM simulations. A few numerical modeling 
methods are observed in literature for modeling PCMs including the enthalpy method, the heat capacity 
method, the temperature transforming model, and the heat source method. The study compares and 
highlights the advantages, disadvantages and limitations of these models and methods. It particularly 
explores the viability of these methods for building applications. The paper further reviews the PCM 
models that have been integrated into prevalent whole building simulation programs such as 
EnergyPlus, TRNSYS, ESP-r etc. The study reveals that the heat capacity method is mostly used in 
programs, despite of its limitations on time and spatial resolutions. Further research is found necessary 
to identify the efficiency and accuracy of these methods in building applications. 
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1. Introduction 


Thermal energy storage (TES) or thermal mass is a property of 
materials that describes its ability to absorb, store and release heat 
depending on the surrounding environmental conditions. Traditional 
architecture, for example, is distinguished with its heavy weight and 
thermally massive construction elements to moderate the indoor 
environment extremes experienced in hot or cold days. The thermal 
properties of construction elements have significantly improved 
thermal comfort by manipulating the indoor air temperature without 
the need of mechanical air conditioning systems [1]. On the other 
hand, light weight buildings are characterised by its lower thermal 
mass and thus expose to significant temperature swings, demanding 
high cooling and heating energy. A dynamic thermal mass such as 
phase change materials (PCM) has been considered as a promising 
technology to reduce the inherited climatic deficiency in light weight 
buildings. 

The apparent advantage of using PCMs lies on the amount of 
latent heat a thin PCM layer can store compared to that in a 
sensible heat storage material such as concrete. For instance, a 
wall of 25 mm thickness with PCM could store an equivalent 
amount of thermal energy as a 420 mm thick concrete wall [2]. 
As a result, the use of PCMs has recently attracted great attentions 
for improving thermal and energy performance of buildings [3-8]. 
Recent studies that review PCMs in building applications can be 
found in [9-18]. Various challenges are, however, arisen when 
using PCM in buildings, including the variety of materials avail- 
able, material liability and safety, cost and economic feasibility, 
design configurations, integration with other sustainable energy 
technologies, impact on thermal and energy performance. The 
problem can then be considered as an optimization dilemma 
where all these counterparts’ challenges are becoming critical in 
the design process. As a result, computational modeling is often 
used as an effective tool to quantitatively understand and help 
resolving this optimization problem. 

This review paper adds to the already existing studies that 
discuss the mathematical modeling of phase change materials for 
building applications such as those described in [19-21]. The 
objective of this study is to systematically review the general 
modeling theories and techniques for PCM, with an emphasis on 
the specific models used for simulating the thermal and energy 
performance of PCMs embedded in building enclosures. Further- 
more, this paper reviews and summarizes the capabilities, 
limitations and validations of prevalent whole building simula- 
tion programs that have been used for modeling phase change 
materials. 


2. General formulation of phase change problems 


The main feature of phase change problems (i.e., Stefan 
problems) is the moving boundary where the Stefan condition 


must be met. For pure materials there is a clear distinction 

between the solid and liquid phase separated by a sharp moving 

interface and hence melting occurs at isothermal temperature. 

For conduction-dominated heat transfer, the governing equation 

can be written for the solid and liquid phase, respectively, which 

have to be satisfied by the Stefan condition as follows [22]: 
Heat transfer in the solid phase: 


oT; ð oT; 
px ex St = 5 (kx mi) o 
Heat transfer in the liquid phase: 
ôT ô aTi 
pxax T= (n x) (2) 


The Stefan condition that enforces the heat balance at the 
solid-liquid interface is: 


a oT, ð ôT, 
Z (kx m) xn 5 (kx a) xn=pxLxvxn (3) 


Very few analytical solutions are available in the closed form 
for phase change problems and can be found in advanced heat 
transfer books such as those by Crank [23], Alexiades and 
Solomon [24], and Özıs ık [25]. Therefore, approximate numerical 
solutions are usually used to handle this class of problems. The 
numerical methods for addressing these problems have been 
reviewed in literature [26-29] and can be generally divided into: 


1. Fixed grid method (i.e., weak solution): These methods consist 
of fixed space grids where the boundary is tracked by the use 
of an auxiliary function. Different approaches are employed to 
account for latent heat evolution [27,29,30]. This class of 
methods has been widely used and therefore will be the focus 
of this paper. 

2. Deforming grid method or front tracking scheme (i.e., classical 
solution or strong numerical solution): These methods allow 
the grid nodes move along with the moving boundary layer 
and thus the space girds deform as the solution develops. 
Here the interface is explicitly tracked using the Stefan condi- 
tion [24]. 

3. Hybrid method: These methods utilize the features of both 
fixed and deforming grids which uses a fixed background grid 
and employs local front tracking schemes to follow the move- 
ment of the boundary [26]. 


3. Numerical formulation of phase change problems using 
fixed grid methods 


An intuitive approach in solving phase change problems is to 
explicitly follow the moving boundary using the front-tracking 
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Nomenclature 

a matrix coefficient 

c specific heat capacity 

f fluid fraction of PCM 

h enthalpy 

k thermal conductivity 

L latent heat of fusion 

n the unit normal on the phase interface 
S source term 

T temperature 

t time 

v velocity of the interface 
x space distance 


Greek letters 


e half range melting temperature 
p density 


Subscripts and superscripts 


A Apparent 

avg average 

e east node 

eff effective 

l liquid state 
m melting 

n iteration level 
p point node 

S solid state 

w west node 


methods. However, this method needs to make a priori assump- 
tion that the boundary is smooth or monotonic during the period 
[31]. This assumption is not always true and therefore reformu- 
lating phase change problems using the fixed grid techniques 
becomes an obvious alternative [23,30,32,33]. The Stefan condi- 
tion Eq.(3) within the fixed grid method is implicitly treated by 
the reformulated governing equation and hence the position of 
the moving boundary is known when the solution is converged. 
The fixed grid method is simple compared to the others, most 
versatile, convenient, adaptable and easily-programmable [24]. 
The latent heat evolution is accounted for in the governing 
equation by using either enthalpy method [34-38], heat capacity 
method [39-42], temperature transforming model [43-46], heat 
source method [38,47-50], or other methods [27,29,51]. The 
following sections will describe the widely used methods. 


3.1. The enthalpy method 


In the Enthalpy method, the latent and specific heat are com- 
bined into an enthalpy term in the governing equation. The enthalpy 
method was proposed by Eyres [38] to deal with variations of 
thermal properties with respect to temperature. For conduction- 
dominated heat transfer, the governing Eqs. (1)-(3) can be reformu- 
lated into one equation where the latent heat is absorbed into the 
enthalpy term as follows: 


oh ə oT 
Paes. (x a (4) 


To demonstrate this method, a fully implicit control volume 
approximation of Eq. (4) for a typical grid shown in Fig. 1 


Fig. 1. A typical control volume grid. 


leads to the following discretized equation: 
h, *' =h; +a! xTw alt x Tp +ag+! x Tit (5) 


According to Eq. (5), it is clear that the current enthalpy (hy! ) 
is dependent on the current value of temperature a and 
therefore the enthalpy term is nonlinear. The equation cannot be 
solved without using proper numerical techniques to handle this 
nonlinearity. This has to be solved either by nonlinear solvers 
such as the Newton’s method or by linearizing the nonlinear 
terms and utilizing iterative methods as fully explained by 
[24,30,31,34,35,52-56]. If a non-linear solver is selected, an 
auxiliary temperature-enthalpy function is required for Eq. (5) 
and can be written for materials that change phases at specific 
temperature range as follows [22]: 


h 
ee hp < Cs x (Tm— €) 
hp + [E+] x(Tm— €) 


p=) ey 


Cix (Tm- €) < hp < Cs x (Tm+ €)+L 


2xe 
hp—(Cs—C)) xTm—L 
— ea 


hp > C x (Tm+ €)+L 
(6) 


Alexiades and Solomon [24] have outlined numerical schemes 
for solving phase change problems with the enthalpy method 
using both linear and nonlinear approaches. Knoll on the other 
hand reviewed various approaches utilizing nonlinear solvers to 
resolve the Stefan problem [55]. He, in particular, developed an 
algorithm to solve the Stefan problem using the Jacobian-free 
Newton-Krylove method and applied for two scenarios: (1) pure 
materials where melting occurs at isothermal temperature and 
(2) non-isothermal case where phase change occurs at a range of 
melting temperature. 

An alternative approach to solving the discretized Eq. (5) is 
to linearize the nonlinear term, hit (T), using the methods 
explained by Patankar [57]. The discretized nonlinear equation 
becomes linear with one primary dependent variable “Temperature” 
that can be iteratively solved with enthalpy using common linear 
solvers such as direct methods (e.g., Gauss elimination or tri- 
diagonal algorithm) or iterative methods (e.g., Gauss-Seidel 
method). Shamsunder [34], for example, proposed a Gauss-Seidel 
iterative scheme where the solution sweeps from west to east to 
determine the state of phase change and subsequently determine 
the new nodal enthalpy. The nodal temperatures are then deter- 
mined based on the discrete form of the enthalpy-temperature 
relationship. To avoid excessive iterations, the scheme was later 
improved by introducing an over-relaxation parameter that is used 
at nodes where no phase change occurs [58]. The scheme was 
however intended for phase change that occurs at isothermal 
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temperature. An iterative Newton linearization scheme was intro- 
duced by Furzeland [31]. The solution process is the same as that of 
Shamsunder except that the over-relaxation parameter can be 
applied at all nodes. 

Iterative methods such as Gauss-Seidel are inherently slow 
and computationally inefficient. Therefore, fast numerical 
schemes have been introduced to improve the computational 
efficiency [35,42,48]. Pham [42] proposed a method that com- 
bines the features of the enthalpy and heat capacity methods. The 
method consists of two steps: a prediction step followed by a 
correction step as shown in Fig. 2. Based on guessed values, the 
new nodal temperatures are predicted (point (2) on the graph). 
The enthalpy is determined based on the predicted temperature 
values. The predicted temperatures are subsequently corrected to 
be consistent with the enthalpy-temperature curve (point (3) on 
the graph). This temperature correction step is the key of this 
method. This method is later known to be the “Quasi-Enthalpy” 
method [59]. 

Voller pointed out that this method might not conserve energy 
at every time step [22] and a better conservative iterative scheme 
was proposed by Swaminathan and Voller as illustrated in Fig. 3 
[35]. The method iterates the predicted and corrected intermedi- 
ate values until the convergence is achieved. The method has 
been recently investigated as an alternative to overcome the 
limitations of the PCM simulation algorithm implemented in 
ESP-r [60]. 


3.2. The heat capacity method 


The heat capacity term in the governing equation imitates the 
effect of enthalpy (sensible and latent heat) by increasing the heat 
capacity value during the phase changing stage. Two approaches 
are generally used to account for the latent heat liberation: the 
apparent heat capacity [22,27,29] and the effective heat capacity 
[61,62]. Although the two approaches differ in the heat capacity 
approximation, recent literatures, however, use the terminologies 
interchangeably. More details on the effective heat capacity 
concept are explained by Poirier [63]. 

The apparent heat capacity method was introduced by 
Hashemi and Sliepcevich [64] to solve a one-dimensional heat 
transfer with phase change in a mushy region. The conduction- 
dominated one-dimensional heat transfer equation using the 


Guess Values (T,,h,) 


Enthalpy, h 


Corrected Values (T,,h,) 


Temperature, T 


Fig. 2. Corrective non-iterative scheme in the Quasi-Enthalpy method at a node 
during one time step. 


Guess Values (T,,h,) 


Enthalpy, h 


Converged Values (T,,h,) 


Temperature, T 


Fig. 3. Corrective iterative scheme in the Enthalpy method at a node during one 
time step. 


apparent heat capacity can be written as: 


am T ðf aT 
oxen xF=5(k 5) (7) 


The method receives the popularity because the temperature 
is the only prime variable that needs to be solved in the 
discretized form. The key in this approach lies in the heat capacity 
approximation. Two methods are commonly used to approximate 
the apparent heat capacity term in Eq. (7): the analytical/empiri- 
cal relationships and the numerical approximations. 


3.2.1. The analytical/empirical relationships 

The heat capacity of a PCM can be determined from the testing 
data with differential scanning calorimeters (DSC). Manufacturers 
of PCMs normally provide limited data pertained to their products 
such as melting temperature, heat of fusion and heat capacity at 
solid and liquid states. Such minimal data can be used to 
approximate the heat capacity of a PCM using a simple direct 
relationship with an introduction of fictitious melting tempera- 
ture range (2 x e) [22,43]: 


Cs, T<Tm-—e€ (Solid region) 
C=} SFG4 Lk  Tm—-e <T<Tm+e (Mushy region) (8) 
ian T>Tm+e¢ (Liquid region) 


Convergence might be an issue when solving Eq. (7), if the half 
phase change range (e) is set too small or the time step is too 
large. There is a possible risk of missing the latent heat contribu- 
tion in a large time step. Hence, DSC testing results can be used to 
form an empirical expression to approximate the heat capacity. 
Fang [65], for instance, proposed a mathematical expression for 
the heat capacity of paraffin—based PCM obtained from DSC. 
Others have suggested and used alternative forms to approximate 
the heat capacity [66-69]. 


3.2.2. The numerical approximations 

Numerical approximation is an alternative when detailed 
information about PCM’s thermal behaviors is available. Many 
numerical approximations have been proposed in literature 
[41,70-75]. For example, Comini [70] applied a numerical tech- 
nique in the finite element method where the heat capacity was 
determined using a derivative of enthalpy with respect to tem- 
perature. Later, Morgan [71] has improved the relationship to 
avoid the convergence problems. When using an iterative scheme, 
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the heat capacity can be approximated using the successive 

temperature and enthalpy solutions. The temporal averaging 

proposed by Morgan [71] is illustrated in Fig. 4 and is represented 

by the following equation: 
_ Ah mth 


g a 


AT T” = Tul (9) 


On the other hand, Lemmon [73] proposed an approximation 
based on the space average rather than the time average 
approach. The temporal and space average approximations are, 
however, prone to convergence issues unless some precautions 
are taken [76]. Solutions to the limitations of the apparent heat 
capacity method have been proposed in literature [33,61,77-79]. 
Voller [22] found that the apparent heat capacity approximation 
based on the direct relationships are more accurate than the 
Morgan approximation used for the cases he studied. 


3.3. The temperature transforming model 


The temperature transforming model was developed by Cao 
and Faghri [80] to overcome the time and spatial limitations in 
the heat capacity method. The model has been used by Faghri and 
his co-workers for many applications [44,45,81]. The method is 
also called “the improved temperature-based equivalent heat 
capacity method” [82]. While the method was tested against 
several benchmark examples, it has been reported to produce 
inconsistent results especially when mass transfer through PCM is 
considered. Corrections were proposed to improve the accuracy 
[81,82]. The key of this method is that the energy Eq. (4) is 
transformed into a nonlinear Eq. (10) with a single dependent 
variable “Temperature” [43]. 
oT ð oT os 

(a) 


p x Cer(T) x a ae pxa (10) 


X 


OX 


This source term is represented by the following Equation [43]. 


Cs; x €, T<Tm- € 
S(T) = (S55) xe + i, Tree <T<Trte (11) 
Cxe +L, T>Tm+ € 


The latent heat during the phase change stage is represented 
by a source term in the governing equation with the heat capacity 
term similar to the apparent heat capacity method. The method 


Guess Values (T,,h,) 


Enthalpy, h 


Converged Values (T,,h,) 


Temperature, T 


Fig. 4. Apparent heat capacity approximation at a node during one time step using 
iterative methods. 


is, however, not commonly used but offers an alternative solution 
when compared to the apparent heat capacity method. 


3.4. The heat source method 


Using the heat source method, the total enthalpy in the 
governing Eq. (4) is split into the specific heat and latent heat 
where the latent heat acts as a source term [23,47]. Eq. (4) thus 
becomes: 


aT af a af, 
px Coe x F = ok T) pxLx (12) 


The method was alluded by Eyres [38] in the mid-1940s. 
In popular schemes, the phase change front is tracked by the 
evaluation of a nodal liquid fraction field which takes a value of 
0 for solid, 1 for liquid, and a value in the range of [0-1] for the 
mushy region [23,49]. With this approach, the fluid fraction is 
linearized and the equation can be solved iteratively with tem- 
perature. The liquid fraction can be approximated using the 
following auxiliary equation [49]: 


0, if T<Tm—e 
fi= STB. if fm—e <T<Tmte (13) 
1, ifT>Tmt+e 


When discretizing Eq. (12) with a fully implicit scheme and 
linearizing the source term “liquid fraction” at the current time 
step, the discretized equation becomes linear and needs to be 
solved for temperature in an iterative manner with the liquid 
fraction. Costa [50] has used this method to numerically simulate 
the latent heat thermal storage. 


3.5. Summary 


Different mathematical models and methods have been sug- 
gested in literature to deal with phase change problems using the 
fixed grids method: enthalpy, heat capacity, temperature trans- 
forming method, and heat source method. Every method has its 
main distinct feature for the latent heat liberation with advantages 
and disadvantages. Table 1 summarizes these methods, and high- 
lights the main feature and their advantages and disadvantages. For 
many reasons including computational efficiency, modeling accu- 
racy and flexibility in selecting solution schemes, the enthalpy 
method is merited to be an attractive mathematical model over 
others for simulating phase change problems. In particular, it 
becomes appealing when the corrective iterative scheme (i.e. a 
fast and energy conservative approach), or non-iterative 
scheme (i.e., a quick but conservative approach at low time steps) 
are implemented. To further exploit these two features for large 
time steps, a quick but energy conservative approach is envisioned. 


4. Models for building enclosures with PCM 


A few models have been developed to solve phase change 
problems on the basis of the general mathematical methods 
described above. A list of the models for various engineering fields 
including building applications was reviewed recently by these 
studies [19,20,83-86]. This section, however, provides a more con- 
centrated and in-depth review on the models that are proposed and 
used for simulating PCM integrated within building enclosures. 

A few innovative and sustainable designs have been proposed 
by integrating PCM within building construction elements. These 
designs demand different levels of model complexity to evaluate 
the thermal performance of such elements. The computational 
models are classified hereafter as the simplified, intermediate and 
sophisticated models. Within this context, the simplified models 
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Feature, advantages and disadvantages of mathematical methods used for phase change problems. 


Mathematical Main feature Advantages Disadvantages Possible solution schemes References 
model for 
latent heat 
evolution 
Enthalpy Enthalpy accounts for e Fast if proper e Difficult to handle supercooling Iterative scheme with Non-linear [24,55] 
method sensible and latent heat scheme is selected problems solvers (e.g., Newton’s methods) 
e Deal with sharp as @ The temperature at a typical grid —_ Linearized-Enthalpy: Corrective [22,35] 
well as gradual point may oscillate with time iterative scheme 
phase change Quasi-Enthalpy: Non-iterative [42] 


Heat capacity Heat capacity accounts Intuitive since 


Lack of computational efficiency 


Temperature correction scheme 


Iterative Scheme (e.g., Gauss- [22,27,29,61,63,64] 


method for both sensible and dealing with one e Small time step and fine grids are Seidel iterative scheme ) if a 
latent heat dependent required for accuracy proper heat capacity 
variable e Difficult in handling cases where approximation is selected 
“Temperature” the phase-change temperature 
e Easy to program range is small 
e Suitable for e Difficult to obtain convergence with 
gradual this technique, and there is always 
phase change a chance that the latent heat is 
underestimated 
e Not applicable for cases where 
phase change occurs at fixed 
temperature 
Temperature Heat capacity and source @ Deal with sharp e Not a common method and Iterative Scheme (e.g., Gauss- [43-46,80,81] 
transforming term are used to account and gradual therefore not tested to evaluate the Seidel iterative scheme ) after 
method for sensible and latent phase change pros and cons linearizing the source term 
heat e Handle large time 
step and 
course grids 
Heat source Latent heat is treated asa œ Intuitive due to e Requires under-relaxation and Iterative Scheme (e.g., Gauss- [23,38,47,49] 


method source term separating the 
latent heat from 
sensible 

Deal with sharp 
and gradual 


phase change 


therefore extra efforts is needed to 
determine the optimum 

relaxation factor 

Lack of computational efficiency 
Problems with round off errors if 
melting occurs over 


Seidel iterative scheme ) after 
linearizing the source term 


temperature range 


are rough approximations of the physics in the phase change 
process but offers quick results. The intermediate models are a 
tradeoff between the speed of the simplified models and the 
accuracy and flexibility of the sophisticated models. The sophis- 
ticated models are created using well validated numerical 
packages that offer a choice of established and optimized numer- 
ical methods. This class offers a high level of accuracy and 
modeling flexibility but is computationally expensive. 


4.1. The simplified models 


Detailed models for simulating PCM within building enclo- 
sures may capture more physics of heat transfer process. How- 
ever, simplified models are sometimes preferred to provide a 
quick estimation of PCM’s thermal performance. Some simplified 
models have been developed with this intention [87-91]. 

A steady-state analytical model for evaluating the benefits of PCM 
in walls and roofs has been proposed by Kaushik [87]. The model 
used the heat capacity method to represent the dynamic thermal 
storage of PCM. The model was utilized to analyze the dynamic 
thermal performance of a free floating building with PCM embedded 
in a south wall façade [88]. The result for a typical mild winter day in 
New-Delhi showed that the wall with PCM outperformed that of an 
ordinary wall. A rough model utilizing the heat capacity method was 
developed to characterize the heat transfer process and subsequently 
estimate the temperature trend in a PCM mixed with gypsum plaster 
board [89]. The simplified model was able to capture the overall 


trend of air temperature in the conditioned room. Another simplified 
physical model using the R-C network method was developed and 
validated for three wall types: light, medium and heavy with shaped- 
stabilized phase change material [90]. The model, however, had 
to use a genetic algorithm to identify the key model parameters: 
resistances and capacitances of the wall layers to reach an optimal 
node distribution. When the optimal parameters are identified, the 
model can be used to simulate the heat transfer process in a wall 
unit that has a PCM layer. Although the model is intended to be 
simple, multiple procedures are necessary for practical applications. 
The model was however implemented to evaluate the energy 
performance of an office building with shaped-stabilized phase 
change material embedded in a wall unit [91]. 


4.2. The intermediate models 


A variety of intermediate models using, respectively, the 
enthalpy method, heat capacity method and heat source methods 
have been developed for one, two and three dimensional cases for 
building enclosures. 


4.2.1. The enthalpy method 

In the enthalpy method, the enthalpy may be solved by 
nonlinear solvers with an auxiliary function (e.g., temperature- 
enthalpy relationship) or implicitly in the governing equation 
using linearization techniques. A theoretical analysis based on the 
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enthalpy method was presented in a study evaluating a PCM in 
a wallboard for solar energy storage [92]. A semi-implicit Crank- 
Nicolson method was used for numerical discretization, which 
was subsequently solved using the Newton’s method. A more 
sophisticated two dimensional finite volume heat transfer model 
based on the enthalpy method was developed and validated to 
explore the behaviors of phase change materials incorporated into 
building-integrated photovoltaic (BIPV) system [93,94]. The heat 
equation is solved with an auxiliary temperature-enthalpy func- 
tion. The model was utilized to perform an optimization study of 
commercially available PCM products embedded into cavity-wall 
systems with different wall-PCM configurations [95]. In addition 
to simulating heat transfer process, the model has the capability 
to solve the Navier-Stokes equation (i.e., the momentum and 
mass equations). The model was expanded later to evaluate a 
three dimensional heat transfer process with PCM [96]. It was 
found that the 3D model does not offer additional accuracy when 
compared to the previously validated 2D model. Another example 
of a validated finite element 3D numerical model based on the 
enthalpy method has been suggested to simulate the PCM mixed 
with common mortars for wall plaster [97]. 

Using the enthalpy linearization approach, a model was 
recently presented as an alternative method for a PCM algorithm 
in ESP-r, a whole building simulation program [60,98]. The 
MATLAB simulation environment was used to develop a one 
dimensional numerical model using a corrective iterative scheme 
proposed by Swaminathan and Voller [35] based on the enthalpy 
linearization. The customized model in MATLAB uses the finite 
volume method with a Crank—-Nicholson scheme to produce a fair 
comparison to ESP-r. The model has proven to be accurate and 
fast when compared to the ESP-r results for a BESTEST Case 600 
model configured with PCM. 


4.2.2. The heat capacity method 

Phase change materials for building applications such as 
Paraffin melt or freeze over a temperature range compared to 
pure materials where phase change occurs at fixed temperature 
[15-17,85]. This property makes the heat capacity method an 
attractive approach to simulating PCM in building applications. 
Utilizing MATLAB package, a research group has developed an 
implicit one dimensional finite difference model for PCM in inner 
wallboard, ceiling and floor with the heat capacity method [99]. 
The discretized equation was solved using the Gauss-Seidel 
iterative method. Although the lab experiments were limited 
and simulation program was incomplete at that stage, the overall 
benefits from PCM in wallboard were evident. 

A semi-implicit one dimensional finite volume heat transfer 
model for simulating PCM in a ceiling of a room using the heat 
capacity method was developed and validated by Pasupathy 
[100,101]. The model was solved using the tri-diagonal matrix 
algorithm (TDMA) with very small time step. Although the overall 
trend of indoor air temperature was captured by the model, the 
numerical results were not in a good agreement with the experi- 
ments due to many limitations of the model. The same model was 
later used for evaluating the PCM integrated into a roof system [102]. 

The heat capacity method has also been implemented in a one 
dimensional numerical model to evaluate shape-stabilized phase 
change materials embedded with a floor heating system [103]. 
The specific heat capacity was used to account for the enthalpy 
of PCM at different temperature regimes. The model gave good 
agreement results when compared to experimental data. The 
model was also used for PCM evaluations under different climates 
and various system configurations [5,104,105]. A variety of 
modeling applications of PCM embedded in floor system for 


different purposes using the heat capacity method have been 
reported in literature [106-108]. 

PCMs have also been integrated in transparent building 
envelopes such as glazed windows. An explicit one dimensional 
finite-difference model based on the heat capacity method was 
extended to evaluate PCM performance when integrated into a 
double glazing system [109,110]. The developed model was 
validated with experimental data and then subsequently utilized 
to evaluate the impact of PCM on heat loss and gains. 

While models developed for building envelope are normally for 
one dimensional geometry, two and three dimensional heat 
transfer approaches have also been suggested for modeling PCM 
using the heat capacity method. In the early 1990s, a numerical 
code “WALL88” was proposed for modeling two dimensional 
transient thermal transport and storage of both sensible and latent 
heat [111]. The model was validated against analytical solution 
and experimental lab results. The model was found to give an 
excellent agreement with experimental results only after allowing 
the PCM to melt over a temperature range rather than at 
isothermal temperature. A three dimensional finite-difference 
heat transfer model using the heat capacity method was devel- 
oped to study the thermal performance of randomly mixed PCM 
and laminated PCM-wallboard systems [112,113]. Although the 
numerical model was not validated in these papers, the simulation 
results were helpful to conclude that laminated PCM-wallboard 
performs thermally better than the randomly PCM-wallboard. This 
model was later validated against an experiment and found a 
maximum of 3% deviation from the average experimental results 
[114]. The model was further used to evaluate the PCM applica- 
tions in the drywall in a passive solar building [115]. The results 
confirmed the conclusions from previous studies for the applica- 
tion of laminated PCM-wallboard. Optimizing the PCM distribu- 
tion within building envelope is the overall goal of an energy 
efficient design. The heat capacity method was adopted by an in- 
house software “CODYMUR”, developed by a team from France, to 
optimize the use of a PCM wallboard for building energy use [116]. 


4.2.3. The heat source method 

The heat source method is an intuitive approach due to the 
separation of specific and latent heat. An explicit one dimensional 
finite-difference heat transfer model for a wall with PCM was 
developed using the heat source method by Athienitis [117] and 
validated against a full-scale outdoor test-room with PCM gyp- 
sum board at interior side. The model showed a reasonable 
agreement with experimental results. A heat transfer model of a 
newly developed hybrid thermal energy storage system (HTESS) 
using PCM capsules in a wall-unit was developed and validated 
for managing solar and electric energy [118]. The numerical 
model uses the heat source method similar to that proposed by 
Voller [48] where the latent heat evolution is represented by a 
source term in the governing equation. The fluid fraction is the 
key to track the latent heat process. 

Phase change materials incorporated within floor systems was 
evaluated using a one dimensional finite volume heat transfer 
model based on the heat source method [119]. The numerical 
model was validated with a benchmark analytical solution of 
Stefan problem explained by Hu [29]. After optimizing the grid 
and the time resolution, the developed model together with an 
optimization algorithm was used to perform an optimization 
analysis on PCM-floor designs. A two dimensional numerical 
model was developed based on the heat source method to 
simulate the effect of PCM in the design of a solar passive wall 
[120]. The model has been verified using benchmark cases 
documented in literature. The model was later validated using 
experimental data performed in the Lab and found to be 
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unsatisfactory due to the limitation of handling the super-cooling 
effect inherited in the tested PCM [121]. 


4.3. The sophisticated models 


Developing a numerical model in two or three dimensional domain 
is complicated and difficult to be generalized for different geometries, 
applications and physics; hence existing simulation packages such as 
COMSOL (formerly known as FEMLAB) [122], ANSYS-FLUENT [123], 
HEATING [124] and others are used as convenient design tools. 
Although, these models offer high level of flexibility, they are not 
fully explored for heat transfer process with phase change. 

One study has used a commercial package FEMLAB (later 
COMSOL multiphysics) to develop a wall model with phase 
change materials using the enthalpy method and heat capacity 
method [125]. COMSOL is a finite element simulation package 
that allows multi-physics modeling for many engineering appli- 
cations. Utilizing this package, the created model was validated 
against experimental results. It was found that both numerical 
methods give good estimation of latent heat evolution process. 
However, the heat capacity method found to be more precise with 
the experimental results when a narrow melting temperature 
range of 2 °C was selected. COMSOL has also been used to study 
envelope systems with PCM [126]. The numerical results from 
COMSOL were successfully compared with another well- 
established numerical model “WUFI-5”. COMSOL is flexible in 
modeling multi-physics within irregular and complex geometries. 
For example, an innovative honeycomb wallboards with PCM 
have been modeled in a 3D domain using the heat capacity 
method [127]. The simulation results showed a very good agree- 
ment with the experimental results. 


Table 2 
Modeling approaches for latent heat evolution in building enclosure. 
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A computational fluid dynamics (CFD) simulation package 
“FLUENT” has been utilized to evaluate a heat source method 
when using user-defined functions for heating and cooling cycles 
of PCM rather than using one idealized function to represent both 
phenomena [128]. The results of a PCM box model utilizing the 
two functions showed very small error (in root mean square 
(RMS) values) when compared with a case of using an ideal 
function for phase change. Another example of using FLUENT is 
reported for PCM integration into a wall cavity system [129]. 

Heat Engineering And Transfer In Nine Geometries (HEATING) 
is a multidimensional, general-purpose heat transfer code that 
has been extensively validated under ASHRAE project RP-1145 
[130]. The code can also be used to model the phase change using 
the heat capacity method. Ahmad [131] has used the program to 
study the behaviors of PCM in wallboard of a test cell. The model 
was validated using experimental test results and found to agree 
well with experiments. The PCM research program at Oak Ridge 
National Laboratory (ORNL) has used this program to study the 
thermal behaviors of PCM in complex two and three geometries 
in building envelope [132]. Lab tests using a heat flow meter 
apparatus (HFMA) have been conducted to validate the model in 
HEATING. HEATING has also been used as a standard benchmark 
numerical package to validate the finite-difference algorithm 
used for PCM modeling in EnergyPlus [133]. 


4.4. Summary 


A variety of models for different building enclosures have been 
developed using various simple, intermediate and sophisticated 
approaches. Table 2 summarizes the models, application usages, 
and validations. It is obvious that very few simplified models have 


Complexity Latent heat evolution’s approach Building enclosure Modeling Solution strategy Validation References 
level case studied formulation 
Simplified models 
Heat capacity method Wall and roof Steady-state 87,88] 
analytical model 
Wallboard Experimental 89 
Optimum nodes for heat capacity Wall R-C Network N/A 90,91] 
distribution using genetic algorithm 
Intermediate models 
Enthalpy method Wall FVM: 1D Newton’s method 92 
BIPV FVM: 2D & 3D Non-linear solver Experimental 93,95,96] 
Wall FVM: 1D Iterative corrective Comparative 60,98] 
scheme 
Heat capacity method Wall FDM: 1D G-S Experimental 99 
Ceiling/roof FVM: 1D TDMA Experimental 100-102] 
Floor FDM: 1D G-S Experimental 5,103-105] 
Glazed-Windows FDM: 1D EM Experimental 109,110] 
Wallboard FDM:2D Analytical and 111 
experimental 
Wallboard FDM: 3D Experimental 112-115] 
Heat source method Wall FDM: 1D EM Experimental 117, 
Wall FDM: 1D TDMA Experimental 118 
Wall FVM: 1D Iterative scheme Analytical 119 
Wall FVM: 2D TDMA Analytical, comparative 120] 
and experimental 
Sophisticated numerical packages 
COMSOL Heat capacity method and heat source Wall FEM: 2D Experimental 125 
method 
Heat capacity method Wall FEM Comparative 126 
Heat capacity method Wall FEM: 3D Experimental 127 
FLUENT Heat source method Wall FVM: 3D SIMPLE algorithm Experimental 128 
Heat source method Wall FVM: 2D Experimental 129 
HEATING Heat capacity method Wall and roof FDM: 1D, 2D and Point-successive over- Experimental 131-133] 


3D 


relaxation iteration 


FVM: finite volume method, FDM: finite difference method, FEM: finite element method, G-S: Gauss-Seidel iterative method, TDMA: tridiagonal matrix algorithm, 


EM: explicit time stepping marching, R-C: resistance-capacitance. 
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been suggested due to the complexity of approximating the heat 
transfer process associated with phase change. The intermediate 
models are commonly used but are developed for specific appli- 
cations and to investigate explicit envelope designs. Hence, they 
lack flexibility in analyzing complex and advanced design alter- 
natives which becomes norms for selecting optimal or near 
optimal designs. Sophisticated models offer flexibility in solving 
complex and multi-physics problems but are not fully explored 
for modeling PCMs. This is partly due to the computational 
inefficiency. They additionally demand considerable amount of 
detailed data inputs, lengthy model setup and validations, and 
limited access to the source codes. 

Generally, all models that adopt the heat capacity or heat 
source methods must be, however, used with small time steps to 
attain acceptable accuracy and therefore slow for whole year 
simulation which is typical for building’s thermal performance 
evaluation. In addition, many existing models ignore inherited 
characteristics of some PCMs such as hysteresis or subcooling and 
therefore cannot be used for this particular application. 


5. PCM models integrated into whole building simulation 
programs 


Many detailed simulation programs are nowadays available to 
assist designers, researchers, manufacturing companies to imple- 
ment new technologies and evaluate innovative ideas that 
improve the energy and thermal performance of buildings. 
Detailed simulation tools perform computations on an hourly or 
sub-hourly bases for accurate considerations of the dynamic 
interactions between all thermal-based elements associated with 
comfort and energy consumption, including building envelope, 
HVAC systems, lighting and control devices [134]. Many building 
simulation tools are listed at the U.S. Department of Energy (DOE) 
web directory [135]. The twenty prevalent whole building energy 
simulation programs that are considered accurate and capable of 
handling the dynamic behaviors of a building and its systems are 
reviewed by Crawley et al. [136]. Few whole building simulation 
programs can handle the thermal performance of building envel- 
ope with phase change materials such as EnergyPlus, TRNSYS, 
ESP-r, and BSim. In addition, some other programs with limited 
capabilities are available for modeling phase change in buildings. 
The following paragraphs brief and compare the conditions of 
these programs. 


5.1. EnergyPlus 


EnergyPlus uses the Conduction Transfer Functions (CTF) to 
approximate heat transfer in building envelope. Since the CTF 
method uses the historical values of heat flux in the computation, 
Barbour [137] has studied the possibility of using this method in 
EnergyPlus to approximate the latent heat evolution in building 
envelope. The study developed multiple sets of CTFs based on the 
temperature of phase change materials. A switching mechanism 
was proposed to exchange between these sets during the simula- 
tion. The CTF-switching algorithm was found to be within 20% 
accuracy for a range of conditions typically encountered in 
buildings. 

The capability of modeling PCMs has been facilitated in 
EnergyPlus program Version 2.0 released in April 2007 by adding 
a conduction finite difference (CondFD) solution algorithm [138]. 
The algorithm uses a semi-implicit finite difference scheme based 
on the heat capacity method with an auxiliary enthalpy- 
temperature dataset to account for latent heat evolution [139]. 
Using this dataset, the heat capacity is approximated using a 
temporal averaging approach similar to that proposed by Morgan 


[71]. While the previous versions of EnergyPlus had a semi- 
implicit scheme for modeling PCMs, a fully implicit scheme has 
been recently added to Version 7 of the program with more 
numerical flexibility [140]. For both schemes, it is however 
recommended to use a small time step for accurate results. 

Experimental validations have been conducted for this algo- 
rithm with mixed feelings of accuracy. Castell, for example [141], 
found that EnergyPlus simulation results did not show a good 
agreement with the experiments when phase change materials 
were implemented in concrete blocks. The study concluded that 
the simulation results did not reflect the thermal improvement of 
PCM observed in the test cells. However, the study highlighted 
that the weather data used in this simulation was not represen- 
tative of the actual weather data. 

An experimental test shed with a commercial PCM product has 
been used to validate EnergyPlus (Version 5) simulation results 
under the climatic conditions of Phoenix, Arizona [142]. It was 
found that the predicted energy consumptions were half during 
winter and slightly greater for the summer months. In addition, 
the time shift was observed for a very short time span during the 
month of April (3 min) and October (9 min). 

Under the climatic conditions of Auckland New-Zealand, an 
experimental study using PCM in gypsum board has been com- 
pared to the simulation results of EnergyPlus using both historical 
weather data and actual measured weather data [143]. Although 
EnergyPlus model using actual weather data has captured the 
overall trend of indoor air temperature but failed to accurately 
predict the actual indoor air temperature from measurements. 
The study highlighted that due to many parameters including air 
infiltration, the simulation results might deviate from the actual 
measurement. 

An early successful validation of the CondFD solution algo- 
rithm used for PCM modeling has been reported by Zhuang [144] 
using two envelope systems with PCM: envelope “A” (one layer of 
PCM: melting temperature at 40 °C) and envelope “B” (two layers 
of PCM: one melting temperature at 40 °C and another at 33 °C). 
The study shows that the largest relative difference in indoor 
temperature is 12.41% and the least is 0.71% between the 
simulation and testing results in a sequential 36h period on 
envelope “A” condition. For envelope “B” condition, the largest 
relative difference is 8.33% and the least is 0.33% in a sequential 
72 h. It was concluded that the most important factors in redu- 
cing the discrepancies between the simulation and the test results 
are to use proper actual weather data as well as using proper 
material thermal characteristics. Other successful validations of 
EnergyPlus algorithm for PCM have been conducted by Campbell 
[145] and Chan [146] using published experimental data by 
Kuznik [147]. For both validation studies, the indoor air tempera- 
ture was found to agree well with the experimental results. 

A field test house was used to study the impacts of PCMs in 
building envelope and consequently used to validate the CondFD 
algorithm in EnergyPlus [148]. The test house used the cellulose 
insulation mixed with 20% PCM by mass. The study reported that 
simulated daily average heat flux through walls was within 9% of 
the field measurements. In addition, simulation results for tem- 
perature distribution through envelope compared fairly well with 
the experimental data apart from some delayed response com- 
pared to the measurement. However, EnergyPlus has given 
unreasonable results for heat fluxes and temperature distribu- 
tions in the attic floor of the experimental house. 

In addition to the validations above, the EnergyPlus’s devel- 
oper team has performed rigorous validation and verification 
studies for general heat transfer calculations as well as the 
CondFD solution algorithm [133,149,150]. These validation stu- 
dies used analytical benchmark solutions, comparative tests with 
well-established program “HEATING”, and experimental results. 
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The studies concluded that versions prior Version 7 contains two 
bugs and subsequently will be fixed in a later version. It is also 
recommended to follow guidelines and bear in mind limitations 
in using Version 7: 


e The time step should be shorter than 3 min. 

e For an accurate hourly thermal performance, 1/3 of the default 
node space should be used. 

e Hysteresis in PCM is not modeled in EnergyPlus and therefore 
inaccurate results may be produced. 


5.2. TRNSYS 


TRNSYS is a modular program where components modules 
“TYPES” are linked together in which output of one type can be 
an input to another in the model. It has been widely used for 
modeling building and its complex systems. Due to its modular- 
ity, users can either utilize available types in the simulation 
package or develop new modules and easily integrate to the 
TRNSYS simulation package. Many features have been introduced 
in Version 16, including a graphical user interface “Simulation 
studio” and the possibility to call external programs such as 
MATLAB, FLUENT and many others [151]. 

Many models have been proposed in TRNSYS for modeling 
phase change heat transfer in building envelope but majority are 
proprietary research modules. Ghoneim for example used a 
modified type of the thermal storage wall “TYPE36” where the 
use of PCM has been tested for thermal storage in a wall system 
[152,153]. The model was based on the enthalpy method and 
solved using an explicit scheme. Despite the numerical problems 
encountered when modeling PCM due to the smaller time step 
required for the stability of the explicit scheme, the model was 
successfully integrated into TRNSYS and validated against pub- 
lished data for a concrete storage wall. Another explicit numerical 
scheme using the enthalpy method was developed for modeling 
the effects of integrating PCM into a solar wall [154,155]. The 
module “TYPE58” was integrated into TRNSYS to explore 
the significance of heating outside air for ventilation purposes 
in the experimental house. 

Modeling PCM in TRNSYS is recently provided through 
“TYPE204” by a team of researchers from Helsinki University of 
Technology [156]. The model simulates heat transfer through a 
wall in a three dimensional domain using the Crank—Nicolson 
scheme with 729 nodes. The model can indeed use a fully implicit 
or explicit scheme with an appropriate selection of a parameter 
that switches between different schemes. The model uses the 
heat capacity to account for latent heat evolution in the wall. 
Although this type has not been validated in its 3D form due to its 
poor computational efficiency, Ahmad [157] has converted the 3D 
model into a 1D module “TYPE101” and validated the modified 
code. The simulation results were compared to experimental 
results from two test cells: one without PCM and another with 
PCM. While the model without PCM works well when compared 
to the experimental results, the model with PCM overestimates 
the daily peak indoor temperature in the cell. The authors 
outlined several reasons for this discrepancy including: (i) eva- 
luation of the energy transmitted through the window, (ii) 
imprecision in the melting temperature range taken in the heat 
capacity definition, (iii) values of the convective heat transfer 
coefficient between wall surfaces and ambient air and (iv) 
existence of cold bridges. Out of these reasons, it was found that 
correcting cold bridges by introducing extra term for resistance 
improved the simulation results significantly when compared to 
the experimental results. 


A study reported a simplified approach of simulating PCM in 
walls/ceiling and floor in TRNSYS [158]. The approach is to use the 
existing capability of TRSNYS to simulate a standard active wall in 
“TYPE56” (i.e., building module in TRNSYS). The key in this 
approach is a user input of equivalent heat transfer coefficients 
introduced in each time step of the simulation that characterizes 
the thermal behaviors of the wall with PCM. The model does not 
evaluate the real heat transfer behaviors in PCM but accurate 
enough for modeling PCM thermal behaviors. The model has been 
validated under laboratory setting conditions. 

On the other hand, Schranzhofer et al. [159] have developed a 
PCM module “TYPE241” where PCM was modeled as an internal 
layer based on the heat source method. TRNSYS capability was 
utilized to model other envelope layers using the transfer func- 
tion method by creating dummy contact zones between the PCM 
layer and the remaining layers. In this type, the PCM is modeled 
using external code based on finite different method with other 
layers modeled through CTF algorithm available internally in 
TRNSYS “TYPE56”. One advantage of this approach is the short 
computational time needed for numerical solution but the phy- 
sics might not be captured well because of assumptions involved 
in the dummy contact zones. The model however was not 
validated due to a lack of appropriate experimental data. 

Kuznik et al. have recently developed a new model “TYPE260’ in 
TRNSYS utilizing the heat capacity method [160]. The model is semi- 
implicit since the physical properties of PCM used in the computa- 
tions are calculated from previous time step. This type has been 
validated with two lab tests conducted by authors: one when the 
outdoor temperature was increased in two steps and the second 
when it was a sinusoidal behavior. The heating heat capacity curve 
was used for the numerical modeling. For both validations, the 
simulation results showed good agreement with the test results. 

A newly developed one-dimensional heat transfer model using 
the heat capacity method was applied to a dividing wall with 16 
glass bricks filled with PCM in TRNSYS [161]. The model has been 
validated and showed fair agreement with experimental results. In 
addition, a simplified PCM module “TYPE1270” has been recently 
developed by Thermal Energy System Specialists (TESS) and added 
to its commercially available individual components [162]. The 
module simulates PCM as an internal layer within an envelope 
system. The model is currently limited to materials that melt/freeze 
at isothermal temperature and with constant specific heat at solid 
and liquid. In the transition state, the PCM layer temperature is 
constant and the model tracks the energy absorption and release. 
The tracking methodology is similar to the heat source method and 
therefore can be identified as “Quasi-Heat Source Method”. 


5.3. ESP-r 


ESP-r is a dynamic energy simulation tool of UK, used for 
modeling thermal, visual and acoustic performance of buildings 
[163]. With many features suitable to model advance sustainable 
energy technologies, ESP-r has the capability to model phase 
change materials using two methods: the effective heat capacity 
method and the additional heat source method [62,164,165]. ESP- 
r uses four models for PCM simulation, with one that accounts for 
sub-cooling, using special materials function. However, it is 
necessary to use a small time step to obtain accurate results for 
these two methods. While simulation results using ESP-r have 
been found in literature, none showed any substantial validations 
for these two algorithms in ESP-r [166-169]. 


5.4. BSim 


BSim is a dynamic simulation program originated from 
Denmark that offers an easy user graphical interface [170]. Using 
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the quasi-steady state in building modeling, the program models 
phase change using the heat capacity method [171]. The simula- 
tion time step has to be small, too, for accurate prediction. Lab 
test results from literature were used to validate the model on 
three cases: continuous heating, continuous cooling, and heating 
but with initial temperature below melting point of PCM. The 
simulation model captures the overall trend of actual thermal 
behaviors of PCM but with small deviations. 


Table 3 
Numerical methods for latent heat evolution in building simulation programs. 


5.5. Other building simulation programs 


Some other building simulation programs have been devel- 
oped for specific research purposes of modeling phase change 
materials, such as RADCOOL [172], ESim [173], and CoDyBa [66]. 
Few have been proposed and developed to simulate simple 
building configurations, such as PCMExpress [174] or the one 
using Engineering Equation Solver (EES) [106]. Due to limited 


Building Module Numerical Numerical Time stepping Limitations/Constrains Validation Reference 
simulation identification formulation method used for scheme 
latent heat 
evolution 
EnergyPlus CondFD FDM: 1D Heat capacity 1. Implicit e Time step <3 min Analytical, [133,144-146,148-150] 
method 2. Semi-implicit e Small grids Comparative and 
e Hysteresis in PCM is not Experimental 
modeled 
e Phase Change at isothermal 
temperature is not modeled 
TRNSYS Modified FDM: 1D Enthalpy method Explicit e Low time step Limited validation [152,153] 
“TYPE36” e No access to the code using 
experimental 
results for 
concrete 
“TYPE58” FDM: 2D Enthalpy method Explicit No access to the code Experimental [154] 
“TYPE204” FDM: 3D Heat capacity Select an appropriate Computationally inefficient N/A [156] 
method factor for implicit, 
semi-implicit or 
explicit 
“TYPE101” FDM: 1D Heat Capacity Semi-implicit e A correction factor to Experimental [157] 
Method (Crank-Nicolson) account for cold bridges has 
to be used for model 
accuracy 
TRNSYS Equivalent Variable heat e Real heat transfer physics in Experimental [158] 
“Active Wall” heat transfer source function PCM is not modeled 
coefficients mimicking PCM 
behavior 
“TYPE241” FDM: 1D Heat source No Published data N/A [159] 
method 
“TYPE260” FDM: 1D Heat capacity Implicit e Thermal properties including Experimental [160] 
method heat capacity are based on 
previous time step (i.e., 
explicit scheme) 
Modified FDM: 1D Heat capacity Implicit e Developed for Internal Experimental [161] 
“TYPE101” method partition wall 
“TYPE1270” Lumped Quasi-heat source e Very simplified model N/A [162] 
method method e Internal layer within an 
using heat envelope 
balance e Based on lumped heat 
balance (not a finite volume), 
low accuracy 
e For phase change at fixed 
temperature 
ESP-r SPMCMP53- FDM: 1D Heat capacity and e Low time step N/A [164,165] 
SPMCMP56 heat source 
method 
BSim FVM: 1D Heat capacity Implicit e Low time step to avoid Experimental [171] 
method instability 
RADCOOL FDM: 1D Heat capacity Implicit [172] 
method 
ESim FDM: 1D Heat capacity Explicit e Explicit scheme requires low Experimental [179] 
method time step to avoid instability 


FVM: finite volume method, FDM: finite difference method. 
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literature available for these programs, only few will be discussed 
hereafter. 


5.5.1. RADCOOL 

RADCOOL is a design tool for cooling and heating system 
developed at the Lawrence Berkeley National Laboratory of the US 
[175]. The program was created using the simulation problem 
analysis and research kernel (SPARK) [176]. A one dimensional 
finite-difference model for a wall with PCM was added to this 
validated thermal building simulation program [172]. The model 
was then used to study the capability of a double PCM-wallboard 
to achieve thermal comfort without using mechanical cooling 
system under a typical climatic condition of Sunnyvale, California 
[177]. 


5.5.2. ESim 

ESim was developed at University of Dayton for building 
energy simulation and can be downloaded from the developer 
website [178]. The simulation program was expanded and vali- 
dated to model PCM-wallboards using an explicit finite-difference 
approach [179]. A list of template files are available for use but 
with limited capability to model complex buildings and systems. 


5.5.3. PCMExpress 

PCMExpress is a planning and simulation program for build- 
ings using phase change materials [174], which was developed by 
a German company, in collaboration with the Fraunhofer Institute 
for Solar Energy (ISE) in Freiburg and partners from industry 
[180]. The program simulates a free floating building with a 
library data for weather and various construction materials 
including PCM and the flexibility to add new materials. It is an 
effective tool to evaluate the economic and technical feasibility of 
PCM usage during an early design stage. The mathematical model 
of the heat transfer process in PCM is not available. The model has 
been tested by Castell [141] who found that the simulations 
deviate significantly from the experiments. As commented by 
Castell, the discrepancy could be attributed to the lack of accurate 
infiltration model in the program. The program however has been 
used to demonstrate the impact of using PCM in residential and 
commercial buildings [180,181]. 


5.6. Summary 


Whole building simulation programs play an important role 
for studying the economic and technical feasibility of PCMs. This 
section reviews the capability of various whole building simula- 
tion programs as summarized in Table 3. It is noted that most 
PCM models integrated into whole building simulation programs 
are based on the heat capacity method. Hence, it becomes 
necessary to reduce the typical one hour time step to a very 
small time step (i.e., in order of minutes) to achieve acceptable 
level of accuracy. For one year thermal performance evaluation, 
building simulation programs become computationally inefficient 
since iterative methods are used in each time step. Additionally, 
the convergence may not be achieved due to numerical instability 
especially when PCM enters or leaves the phase change region. 
With all constraints and limitations above, none of whole building 
simulation programs are currently implementing efficient math- 
ematical models that are quick, accurate and numerically stable 
at realistic time step. It becomes important to thoroughly inves- 
tigate different mathematical models with various numerical 
approaches for modeling PCMs in whole building simulation 
programs. 


6. Conclusions 


Significant heat storage offered by phase change materials is 
promising and favorable for developing various innovative building 
energy technologies. To quantify the technical and economic feasi- 
bility of PCM-embedded technologies, it requires the development 
of proper computational models. This study reviews various numer- 
ical modeling approaches of phase change problems such as the 
enthalpy, heat capacity, temperature-transforming and heat source 
methods. The main features, advantages and disadvantages of each 
method have been discussed. The discretized form of the heat 
equation with PCM can either be solved with nonlinear solvers such 
as Newton’s methods or via linearizing the nonlinear term and using 
linear solvers such as iterative methods. For both approaches, the 
numerical solutions are computationally inefficient or difficult to 
reach convergence. Therefore, fast numerical schemes are suggested 
such as the quasi enthalpy non-iterative scheme or the enthalpy 
conservative iterative scheme. 

Using these general mathematical methods, different compu- 
tational models have been developed to simulate PCMs in build- 
ing enclosures. Based on the level of complexity, models are 
classified into three categories: simple, intermediate and sophis- 
ticated models. Majority of these models have been validated 
using analytical solutions, comparative testing using validated 
numerical models, and/or experimental results. 

While many models are used to study the heat transfer in an 
enclosure unit, a few models have been integrated into whole 
building simulation programs. A variety of models are available 
and some are available with no cost to users such as EnergyPlus, 
“TYPE204” in TRNSYS or ESP-r. These particular models, however, 
have limitations on modeling PCM including the time and spatial 
resolutions, inability to model hysteresis, lack of validations of 
some models, and poor computational efficiency. These modeling 
challenges add complexity to the already existing uncertainties in 
experimental results of PCM’s thermal behaviors. Therefore, 
further research is needed to quantitatively explore the prediction 
performance of different models including their limitations on 
accuracy, parameters sensitivity, speed, and stability for modeling 
PCM envelopes under different climatic and operating conditions. 
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